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Minimization  Of  Squared  Angular  Error  To  Reduce  Bias 


SUMMARY 

The  systems  reviewed  to  date  use  the  weighted  perpendicular  method.  As 
the  number  of  bearings  increases  the  bias  of  this  approach  becomes  a  larger 
and  larger  percentage  of  the  total  error.  This  does  not  happen  to 
minimization  of  angular  error.  On  the  other  hand,  minimization  of  angular 
error  is  harder  to  compute.  How  much  harder? 

I.  Number  of  Iterations  after  a  Perpendicular  Method  Start- 

Weighted  Perpendicular  converges  in  one  or  two  iterations 
whereas  Minimization  of  Squared  Angular  Error  of  angular  error 
should  probably  use  two  to  three  iterations  (after  using  the 
Perpendicular  Method.) 

II.  Shared  calculations- 

In  the  Weighted  Perpendicular  reweighting  can  be  done  on  saved 
terms  without  significantly  reducing  recalculation.  Minimization 
of  Squared  Angular  Error  must  recalculate  all  terms  each  iteration. 

III.  Number  of  multiplications  and  additions- 

Minimization  of  Squared  Angular  Error  has  slightly  more  than 
Weighted  Perpendicular.  This  difference  should  be  fairly 
insignificant  for  total  program  run  time  however. 

IV.  Computation  of  Transcendental  Functions- 

Minimization  of  Squared  Angular  Error  requires  computing  an 
Inverse  trig.  All  other  methods  require  computation  of  trig 
functions  but  minimization  of  angular  error  does  not.  The  primary 
difference  is  probably  not  the  difference  between  trig  functions 
and  inverse  trig  functions.  The  primary  difference  is  the  trig 
functions  need  only  be  computed  once  whereas  the  inverse  trig 
functions  must  be  recomputed  each  time  a  new  fix  is  found.  To  get 
the  benefits  of  a  fix  with  bias  properties  similar  to  Minimization 
of  Squared  Angular  Error  without  having  to  recompute  transcendental 
functions  one  should  use  Minimization  of  Sine  Squared  of  Angular 
Error.  According  to  its  documentation  this  sin  variation  was  used 
by  FFIX.  MARC  has  not  seen  the  algorithm  FFIX  used.  Our  own  code 
for  this  method  is  included  later  in  this  report. 

MARC  recommends  that  slowdown  In  the  number  of  fixes  that  can  be 
processed  be  compared  against  the  percentage  of  error  resulting  from  fix 
algorithm  bias.  Analysis  of  data  from  field  trials  should  be  an  important 
part  of  determining  this  percentage.  Analysis  based  on  models  can  predict 
what  the  result  will  be  but  only  field  trials  can  validate  that  important 
sources  of  error  have  not  been  left  out  of  a  model.  Models  should  be  a 
consideration  in  the  planning  and  evaluation  of  trials  but  they  are  only  one 


element  of  an  Integrated  investigation. 


The  remainder  of  this  memo  consists  of: 

I.  a  section  giving  a  general  background  on  the  issue  of 

error  from  the  fix  algorithm  (excluding  bearing  selection) 

II.  a  section  giving  the  algorithm  for  one  iteration  of  the 
minimization  of  squared  angular  error 

I.  FIX  ALGORITHM  ERRORS  (EXCLUDING  BEARING  SELECTION  ERROR) 

Analysis  of  fix  algorithms  under  the  standard  normal  angular  error  model 
yields  the  following  conclusions: 

A.  Minimizing  variance: 

Most  of  the  error  for  small  numbers  of  LOBs  is  reflected  in  the 
variance  of  the  estimate.  Minimization  of  the  variance  is  therefore 
critical.  Most  fix  methods  have  variances  that  are  very  similar. 

The  method  of  minimization  of  Perpendicular  Distance  is  an  exception 
and  hence  should  not  be  used.  A  Senior  Thesis  underway  at  Claremont 
McKenna  College  on  a  generalization  of  the  Cramer  Rao  Inequality 
will  show  that  the  variance  of  the  techniques  such  as  Weighted 
Perpendicular ,  Minimization  of  Squared  Angular  Error,  and  Minimization 
of  Sine  Squared  of  Angular  Error  are  essentially  optimal  relative  to 
errors  resulting  from  a  shared  angular  uncertainty.  (If  there  are  known 
differences  in  angular  uncertainty  they  should  be  incorporated.  Any 
of  these  methods  could  be  modified  to  weight  for  proportional  accuracy.) 
The  amount  that  the  variance  of  the  Perpendicular  Method  is  larger 
depends  on  how  much  the  relative  distance  from  the  emitter  to  the 
sensors  varies  depending  on  the  sensor. 

B.  Reducing  Bias: 

For  small  sample  sizes  bias  is  the  second  largest  source  of  error. 

For  large  sample  sizes  it  is  the  largest  source  of  error.  It  can  be 
shown  that  for  most  methods  (certainly  all  Minimization  of  Square  Error 
based  methods)  and  2  bearings  the  different  methods  give  the  same  fix, 
and  hence  have  the  same  bias.  As  sample  size  goes  up  some  methods  have 
bias  shrink  inversely  proportionately  to  sample  size  and  others  do  not. 
The  Perpendicular  Method  and  Weighted  Perpendicular  Method  do  not  have 
shrinking  bias  in  general.  In  fact,  in  some  situations  bias  can 
increase  with  sample  size  for  Perpendicular  and  Weighted  Perpendicular. 

C.  Skewness: 

Intuition  suggests  that  skewness  depends  on  the  relative  size  of  an 
ellipse  compared  to  the  distance  from  the  sensors.  Skewness  should 
decrease  as  sample  size  increases  if  this  is  the  case.  MARC  would 
expect  Weighted  Perpendicular  to  have  skewness  similar  to  Minimization 
of  Squared  Angular  Error.  Furthermore  MARC  expects  that  both  would  be 
somewhat  smaller  than  the  skewness  of  the  unweighted  Perpendicular 
method  in  a  way  directly  related  to  the  larger  variance  of  the 
Perpendicular  method.  MARC  has  not,  however,  tested  these  theories. 

II.  ITERATIONS  TO  UPDATE  LOCATION  ESTIMATES 

Commented  Squared  Angular  Error  Minimization  code  is  Attachment  I. 

Minimization  of  Sin  Squared  of  Angular  Error  code  is  Attachment  II. 

Notes  on  the  code:  X  and  Y  are  used  as  subscripts.  The  initialization 
portion  of  the  program  Is  not  included  but  it  sets  X-1  (for  the  first 
coordinate  or  X  coordinate)  and  Y-2  (for  the  second  coordinate  or  Y 


coordinate.)  ATN  is  the  Arctangent  function.  Bearing(I)  is  the  angle  in 
degrees  from  true  North  of  the  bearing.  The  algorithm  used  is  to  update 
successive  estimates  of  the  fix  by  subtracting  the  inverse  of  the  Hessian 
(second  derivative)  of  the  likelihood  times  the  gradient  of  the  likelihood 
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-  ALGORITHM  LISTING  - 

ITERATION  FROM  MINIMIZATION  OF  ANGULAR  ERROR 


ATTACHMENT  I 


25U8 

25^9 
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2595 
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2597 

2598 

2599 

2600 


!  ********************  Variables  ************************* 

; 

!  Partial_wrt(X)  is  the  partial  of  the  likelihood  L  with  respect  to  X 
!  Partial_wrt( Y)  is  the  partial  of  the  likelihood  L  with  respect  to  Y 
!  L(X,X) ,L(X,Y),  and  L(Y,Y)  are  the  second  partials  of  L 
!  Dif  f*=Angular  Difference 

!  Sensor( I , X)=X  coordinate  of  the  ith  sensor 
!  Sensor ( I , Y) =Y  coordinate  of  the  ith  sensor 
!  Fix(X)=X  coordinate  of  the  current  location  estimate 
!  Fix(Y)=Y  coordinate  of  the  current  location  estimate 
!  Dx  is  the  difference  in  X  coordinates  of  the  estimate  and  the  sensor 
!  Dy  is  the  difference  in  Y  coordinates  of  the  estimate  and  the  sensor 
!  Determinants ( X, X)L(  Y,  Y) ~L( X,  Y) L(  Y , X) 

!  Other  variables  should  be  self-explanatory 

j 

!  ENTER  THIS  LOOP  WITH  AN  ESTIMATE  BASED  ON  THE  PERPENDICULAR  METHOD 

I 

FOR  Iteration- 1  TO  4 

!  ****************  INITIALIZE  TO  COMPUTE  A  SUMMATION  *************** 

Partial_wrt(X)*=0 

Par  t  ial_wrt  ( Y )  -0 

L(X,X)-0 

L(X,Y)-0 

L(Y.Y)  0 

i  *************  COMPUTE  SUMMATION  FOR  L  &  Partial_wrt  ************** 
FOR  I»1  TO  Nsensors 

Dif f-(Bearing(I)-ATN( (Sensor ( I, X)-Fix(X))/ (Sensor ( I, Y)-Fix(Y)))) 

IF  Dif f >180  THEN  Diff-Dif f~360 
IF  Dif f <-180  THEN  Diff-Diff+360 

Diff-Dif f*PI/1 80.0  !  ***  convert  to  radians  *** 

Dx-Fix(X)-Sensor (I, X) 

Dy-Fix( Y)-Sensor(I, Y) 

Dx_2-Dx*Dx 

Dy_2»Dy*Dy 

Dxy-Dx*Dy 

Distance_sq**Dx_2  +  Dy_2 
Distance_!»-Distance_sq*Distance_sq 

Partial_wrt(X)»Partial_wrt(X)-2.0*Dy*Diff/Distance_sq 
Partial_wrt( Y)-Partial_wrt( Y) +2. 0*Dx*Diff/Distance_sq 
L( X, X) -L(X,X)+2.0«( 2. 0*Diff*Dxy+Dy_2) /Distance^ 
L(X,Y)-L(X,Y)+2.0*((2.0*Diff*Dy_2-Dxy)/Distance_l1-Diff/Distance_sq) 
L( Y , Y)-L( Y, Y)+2.0*(-2.0*Dif f *Dxy+Dx_2) /Distance_4 
NEXT  I 

L( Y,X)-L(X, Y)  !  ***  no  sense  computing  it  twice  *** 

j 

!  Now  update  the  Fix  by  subtracting  (L  inverse)  times  (Partial_wrt) 

j 

Determinant»L(XfX)*L(Y, Y)-L(X, Y)*L(Yf X) 

Fix(X)-Fix(X)-(L(Y, Y) *Partial_wrt( X)-L(X, Y) *Partial_wrt( Y) ) /Determinant 
Fix( Y)-Fix( Y)-(-L( Y, X) *Partial_wrt( X) +L( X, X) *Partial_wrt( Y) ) /Determinant 
PRINT  "FIX“(";Fix(X);",";Fix(Y);”)" 

NEXT  Iteration 


ATTACHMENT  II 


-  ALGORITHM  LISTING - 

ITERATION  FROM  MINIMIZATION  OF  SIN  SQUARED  OF  ANGULAR  ERROR 


2629  FDR  1=1  TO  Nsensors 

2630  S2(I)=SIN(2*Bearing(I)) 

2631  C2(I)=CCS(2»Bearing(I)) 

2691  NEXT  I 

2653  FDR  Iterations  TO  3 

2659  !  ****************  INITIALIZE  TO  COMPUTE  A  SUMMATION  a#*#*#**#****** 

2655  Partial_wrt(X)-0 

2656  Partial_wrt(Y)=0 

2657  L(X,X)=0 

2658  L(X,Y)=0 

2659  L(Y,Y)=0 

2662  !  *************  OOMPUTE  SUMMATION  FDR  L  &  Partial_vrt  ************** 

2663  FDR  1=1  TO  Nasnsors 

2666  Dx=Fix(X)-Sensor(I,X) 

2667  Dy=Fix(Y)-Sensor(I,Y) 

2668  Dx_2=Dx*Dx 

2669  Dy_2=Dy*Dy 

2670  Dx_3=Dx_2*Dx 

2671  Dy_3=Dy_2*Dy 

2672  Dx_9=Dx_3*Dx 

2673  Dy_^=Dy_3*Dy 

2679  Dxy=Dx*Dy 

2676  Distance_sq=Dx_2+Dy_2 

2677  D9=Distance_sq*Distance_sq 

2678  D6=Distance_9*Distance_sq 

2679  S=S2(I) 

2680  C=C2(I) 

2681  Part ial_wrt ( X) -Part ial_wrt ( X) +Dy* ( S* ( Dx_2-Dy_2 ) +2 . 0*C*Dxy) /D9 

2682  Par  t  ialjwr  t ( Y ) -Partial_wrt ( Y) +Dx* ( S* ( Dy_2-Dx_2 ) -2 . 0*C*  Dxy ) /D9 

2683  L(X,X)=L(X,X)+2.0*(S*(3.0*Dx*Dy_3-Dx_3*Py)-C*(3.0*Dx_2*Dy_2-Dy_9))/D6 

2689  L(X,Y)=L(X,Y)+(S*(Dx_9-6*Dx_2*Py_2+Dy_9)+9*C*(Dx_3*DrDx«Dy_3))/D6 

2685  L(  Y ,  Y)  =L(  Y ,  Y)  +2 . 0*  ( S*  ( 3  •  0*Dx_3*Dy-Dx*Ety_3 )  +C*  ( 3  •  0*Dx_2*Dy_2-Dx_9 ) )  /D6 

2686  NEXT  I 

2687  L(Y,X)=L(X,Y)  !  ***  no  sense  confuting  it  twioe  *** 

2688  ! 

2689  !  New  update  the  Fix  by  subtracting  (L  inverse)  times  (Partial_wrt) 

2690  ! 

2691  Determinant-L(X,X)*L(Y,Y)-L(X,Y)*L(Y,X) 

2692  Fix(X)=Fix(X)-(L(Y,Y)*Partial_vrt(X)-L(X,Y)«Partial_wrt(Y))/Determinant 

2693  Fix(Y)«Fix(Y)-(-L(Y,X)*Partial  vrt(X)+L(X,X)*Partial_wrt(Y)) /Determinant 

2699  PRINT  "FIX-(";Fix(X);",”;Fix(YT;")" 

2695  NEXT  Iteration 


